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ABSTRACT 

Many optimization techniques have been invented to reduce the noise that is inherent in Monte 
Carlo radiative transfer simulations. As the typical detectors used in Monte Carlo simulations 
do not take into account all the information contained in the impacting photon packages, 
there is still room to optimize this detection process and the corresponding estimate of the 
surface brightness distributions. We want to investigate how all the information contained in 
the distribution of impacting photon packages can be optimally used to decrease the noise 
in the surface brightness distributions and hence to increase the efficiency of Monte Carlo 
radiative transfer simulations. 

We demonstrate that the estimate of the surface brightness distribution in a Monte Carlo 
radiative transfer simulation is similar to the estimate of the density distribution in an SPH 
simulation. Based on this similarity, a recipe is constructed for smart detectors that take full 
advantage of the exact location of the impact of the photon packages. Several types of smart 
detectors, each corresponding to a different smoothing kernel, are presented. We show that 
smart detectors, while preserving the same effective resolution, reduce the noise in the sur- 
face brightness distributions compared to the classical detectors. The most efficient smart 
detector realizes a noise reduction of about 10%, which corresponds to a reduction of the re- 
quired number of photon packages (i.e. a reduction of the simulation run time) of 20%. As 
the practical implementation of the smart detectors is straightforward and the additional com- 
putational cost is completely negligible, we recommend the use of smart detectors in Monte 
Carlo radiative transfer simulations. 
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1 INTRODUCTION 

The Monte Carlo method (e.g. Cashwell & Everett 1959; Bianchi, 
Ferrara, & Giovanardi 1996; Gordon et al. 2001; Baes et al. 2003; 
Niccolini, Woitke, & Lopez 2003) has become one of the most pop- 
ular methods to perform radiative transfer simulations. One of the 
greatest advantages of the Monte Carlo method is its conceptual 
simplicity. Instead of solving the radiative transfer equations that 
describe the radiation field, Monte Carlo simulations actually fol- 
low the photon packages that make up the radiation field in a very 
natural way. This ensures that, compared to other approaches to 
multi-dimensional radiative transfer problems, the practical imple- 
mentation of Monte Carlo radiative transfer is surprisingly straight- 
forward. This simple, straightforward approach also enables the in- 
clusion of additional ingredients, such as the polarization of scat- 
tered light (Code & Whitney 1995; Bianchi, Ferrara, & Giovanardi 
1996), the kinematics of the sinks and sources (Matthews & Wood 
2001; Baes & Dejonghe 2002) and gas ionization and recombina- 
tion processes (Ercolano et al. 2003; Ercolano, Barlow, & Storey 
2005). On the other hand, probably the most important disadvan- 
tage of the Monte Carlo method is the appearance of Poisson noise, 
which is inherently tied to the probabilistic nature of the method. 
In pure Monte Carlo radiative transfer simulations, the noise in any 



observed property goes as l/^/N with N the number of photon 
packages, which means that one needs to quadruple the number of 
photon packages to halve the errors. Due to this slow convergence, 
radiative transfer simulations based on the most simple application 
of the Monte Carlo method are very inefficient and have difficulties 
to compete with other methods. 

Fortunately, in the many years since the first applications of 
the Monte Carlo method to radiative transfer simulations, several 
intelligent optimization techniques have been invented to increase 
the efficiency of the Monte Carlo techniques. These techniques, 
equivalent to so-called variance reduction techniques in Monte 
Carlo integration, aim at reducing the noise by including determin- 
istic elements into the probabilistic simulation. One of the first ex- 
amples of such techniques is the well-known forced first scattering 
technique, already included in the first implementations of Monte 
Carlo radiative transfer (Cashwell & Everett 1959; Mattila 1970). 
Other noise-reduction techniques that have strongly increased the 
efficiency of the Monte Carlo method include the use of weighted 
photon packages (Witt 1977), the peel-off technique (Yusef-Zadeh, 
Morris, & White 1984), the treatment of absorption as a continuous 
rather than a discrete process (Lucy 1999), the frequency distribu- 
tion adjustment technique (Bjorkman & Wood 2001; Baes et al. 
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2005) and the use of polychromatic photon packages (Baes, De- 
jonghe, & Davies 2005; Jonsson 2006). 

One aspect of Monte Carlo radiative transfer where no sig- 
nificant noise reduction techniques have been presented is the last 
step in the life cycle of the photon packages, namely their detec- 
tion. The goal of most Monte Carlo radiative transfer simulations 
is to determine the observed surface brightness distribution at some 
observer's position. To construct the observed surface brightness 
distribution, some kind of detector must be simulated, on which 
the photon packages that leave the system are recorded. In the spirit 
that Monte Carlo radiative transfer simulations mimic the real phys- 
ical processes as naturally as possible, the simulated detectors are 
usually natural (idealized) representations of actual CCD detectors. 
They basically consist of a two-dimensional array of pixels, which 
act as a reservoir for the incoming photon packages. When a pho- 
ton package leaves the system and arrives at the location of the 
observer, the correct bin is determined and the luminosity recorded 
in the bin is increased with the luminosity of the photon package. 
At the end of the simulation, the detector is read out like a CCD 
detector and the surface brightness distribution is constructed. 

While this approach seems the most natural way to simulate 
the detection of photon packages in a Monte Carlo simulation, it 
might not be the most efficient. We must be aware that, although we 
are simulating a real detection as closely as possible, we have more 
information at our disposal than real observers. The maximum in- 
formation that a real observer can obtain (in the academic limit 
of perfect noise-free observations and instruments) when imaging 
with a CCD detector, is the number of photon packages that arrive 
in each of his pixels. As numerical simulators, we have at our dis- 
posal the full information on the precise location of the impact of 
each photon package on the detector. It would be a pity to throw 
away this useful additional information just in order to mimic the 
behaviour of a real CCD detector. 

The goal of this paper is to investigate how this additional in- 
formation can be used to decrease the noise in the estimated sur- 
face brightness distributions and hence to increase the efficiency 
of Monte Carlo radiative transfer simulations. In Section 2 we de- 
scribe the classical way of detecting photon packages as a smooth- 
ing process, similar to the averaging processes encountered in the 
smoothed particle hydrodynamics framework for hydrodynamical 
simulations. We use this similarity between both processes to con- 
struct a new set of smart detectors, which aim at curing two draw- 
backs of the classical detectors. In Section 3 we test the accuracy 
and performance of these smart detectors and compare them with 
the classical detector. The results are summarized and discussed in 
Section 4. 



luminosity Lk carried by photon package to the number of previ- 
ously detected photon packages in this particular bin. At the end of 
the simulation, we determine our estimate h{xij) for the surface 
brightness at the position Xij by summing the contribution of all 
photon packages that have been recorded in that bin and correcting 
for the surface of the bin. 
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where 



N^j{x) = I 



' 1 ifxi - |A < X < + |A 

and^, -|A<y<y, + |A 

else. 
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2.2 The link to SPH 

For an interesting point of view on simulating the detection of pho- 
ton packages and the calculation of the surface brightness distri- 
bution, we now shift to another important technique in computa- 
tional astrophysics, namely smoothed particle hydrodynamics or 
SPH (Gingold & Monaghan 1977; Lucy 1977). SPH is a compu- 
tational technique for hydrodynamical simulations in which the a 
fluid is represented as a finite collection of fluid elements. To esti- 
mate the value of a physical field /(r) at an arbitrary position r, 
we must perform local averages over volumes of nonzero extent 
(i.e. over a finite number of fluid particles). The mean or smoothed 
value of /(r), denoted as f^{r) can be determined through kernel 
estimation. 



jjjw{r-r')f(r')Ar', (5) 



where W(r) is the so-called smoothing kernel, which should be 
normalized according to 



/// 



W(r)dr = 1, 



(6) 



and which typically is strongly peaked about r = 0. In particular, 
if the fluid mass distribution p{r) is represented by a flow of N 
particles with mass mfc, we have 



(7) 



and the estimate ps(r) for the mass distribution is given by 



2 SMART DETECTORS 



2.1 The classical detector 



A typical detector in Monte Carlo simulations is based on a realistic 
CCD detector. It consists of a rectangular two-dimensional array of 
bins of dimension A placed on the plane of the sky. We denote the 
centre positions of these bins as Xij = (xi^yj). 



Xi — Xmin ~h A, 
Vj = Vmm + i A. 



(1) 

(2) 



Ps(^) = '^ink W{r - Vk) 



(8) 



The standard interpretation of this formula is to see the fluid ele- 
ments that make up the representation of the fluid not as point-like 
particles, but as particles with a smoothed smeared out mass distri- 
bution pk{r) = mk W{r — Vk). The total density at any position 
r is then the sum of the contribution of the N particles in the fluid. 

Returning back to our simulated CCD detector, we note that 
we can write expression (3) as 



When the /c'th photon package hits the detector at a certain position 
Xk = {xk,yk),^Q find out in which bin it will end, and we add the 



Is{xij) = ^ Lfc W{xij - Xk), 



(9) 
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with W{x) di function defined as 

"if(x+|A)-if(x-|A) 



W[x) 



H{y+\/\)-H{y-\/\) 



(10) 



with H the Heaviside step function. Comparing equations (9) and 
(8) we immediately see a connection between the computation of 
the (three-dimensional) mass density in fluids in the SPH formalism 
and the computation of the (two-dimensional) surface brightness 
in a Monte Carlo radiative transfer simulation. We can make this 
connection clearer by rewriting equation (9) as 



Xij — x') I{x') &x\ 



with 



I{x) = ^Lfc 5{x - Xk 



(11) 



(12) 



Since this latter expression is nothing but the "true" surface bright- 
ness distribution corresponding to N photon packages hitting the 
detector plane (each photon package results in a Dirac delta func- 
tion), equation (11) is the direct analogue of the SPH basic equa- 
tion (5). The bottom-line is that we can interpret the radiation field 
at the plane of the sky as a fluid, which is represented by a fi- 
nite number of smoothed photon packages. Each smoothed photon 
package corresponds to a smeared out surface brightness distribu- 
tion Ik{x) = LkW{x — Xk)- Formula (9) shows that the total 
observed surface brightness distribution at the positions Xij is the 
sum of the contributions of each of these photon packages. 

2.3 Smart detectors 

We have seen that we can interpret the determination of the surface 
brightness distribution in a Monte Carlo radiative transfer simula- 
tion as a smoothing operation similar to the averaging in SPH hy- 
drodynamical simulations. As in SPH simulations, we can consider 
using a different kernel — in the present case this corresponds to a 
different kind of detector. In principle, there is no restriction on the 
shape of the function W{x) apart from the normalization condition 



W{x)dx = 1, 



(13) 



and the requirement that W{x) isa. centrally peaked function. This 
freedom allows us to construct smart detectors which improve upon 
a number of potential disadvantages of the traditional detector. 

A first drawback of the smoothing kernel (10) is that it is not 
isotropic, meaning that it has a preferential direction. Not all pho- 
ton packages hitting the detector at the same given distance from a 
grid point Xij have the same impact on the surface brightness. For 
example, a photon package hitting the detector at 



Xk 



(14) 
, whereas 



will not contribute at all to the surface brightness h{xij] 
a photon package hitting the detector at 

Xk + A, yj + , (15) 

at the same distance, will fully contribute to the estimate of the 
surface brightness at this grid point. It is obvious that this problem 



is solved when we consider circularly symmetric smoothing kernels 
W{x) = W{R). 

A second drawback of the classical detector is that it does not 
take into account all the information that is contained in the impact- 
ing photon package. For every photon package that falls onto our 
detector, we know the exact location of the impact. The only infor- 
mation that a classical simulated detector uses is the bin in which 
this location falls, without any discrimination of the exact location 
within this bin. As argued in the introduction, it would be a pity to 
not use this information just for the sake of simulating a real CCD 
detector. To solve the second problem, we need to look for kernels 
that give more weight to impacts very close to the grid points than 
to impacts at larger distance, which means that we need a kernel 
W{R) that is a monotonically decreasing function of R. 

These are just minor limitations and still leave a lot of room for 
different smart detectors. We can inspire ourselves on the smooth- 
ing kernels that are often used in SPH simulations. The prototypical 
normalized kernel is a gaussian, in two dimensions. 



W{R) 



TTh^ 



exp 



R^\ 



(16) 



The parameter h in this expression, and in all other kernels we will 
discuss, is called the smoothing length. It gives the width of the area 
over which the smoothing is effective (we will later determine the 
optimal value for this parameter). An important drawback of the 
gaussian kernel is its infinite support, meaning that every photon 
package impacting on the detector has a finite contribution to the 
surface brightness at all grid points. Every photon package hence 
in principle requires a summation over all grid points, at most of 
which it contributes an absolutely negligible contribution. 

These efficiency problems are resolved by introducing 
smoothing kernels with a compact support. Several compact sup- 
port kernels have been used in the SPH literature, the most popular 
of them being kernels based on B- splines. The M3 spline kernel 
(Hockney & Eastwood 1981) in two dimensions takes the form 



W{R) = 



16 



137r/i2 



■ 2u^ 







if < 1^ < I 
if 1 
else 



2 <^< i 



(17) 



with u — R/hAi has a continuous first derivative. A similar spline 
kernel but with a continuous second derivative is the M4 spline 
smoothing kernel, introduced by Monaghan & Lattanzio (1985). 
Its two-dimensional form is 



W{R) 



10 



(2- 



if < 1^ < 1, 
if 1 < n < 2, 
else. 



(18) 



In principle, we have not requested that the kernel is everywhere 
non-negative. Some SPH simulations have adopted so-called su- 
perkernels, which are accurate to third or fourth order, but which 
necessarily become negative in some part of the domain. The most 
famous example of such kernels is the supergaussian kernel (Gin- 
gold & Monaghan 1982) but examples with finite support have also 
been considered (Monaghan 1985, 1992; Capuzzo-Dolcetta & Di 
Lisio 2000). Using such superkernels could potentially increase the 
accuracy of the smoothing process, but it can have serious conse- 
quences when there is a sharp change in the surface brightness: an 
undershoot occurs and the recorded surface brightness may become 
negative. In order to avoid such situations, we stick to positive ker- 
nels. 
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Figure 1. The three smoothing kernels used for the construction of the smart 
detectors. The normalization and the width is the same for all three kernels. 



Figure 4. The total noise parameter IZ, defined in equation (26), of the 
Plummer model for different detectors as a function of the total number of 
photon packages in the simulation. 



2.4 Determination of the smoothing length 

Based on the three smoothing kernels (16), (17) and (18) we can 
construct three different kinds of smart detectors. Remains to de- 
termine which value to take for the smoothing length h for each of 
these smart detectors. Obviously, h should be large enough to make 
the smoothing operation meaningful, whereas too large values of h 
will tend to over- smooth and wash away the details of the surface 
brightness distribution. We want the smoothing lengths comparable 
with the pixel scale, such that the noise in the images remains un- 
corrected on a pixel-by-pixel scale (as for a classical detector). Our 
determination of the optimal smoothing length is derived from the 
requirement that the effective area or "resolution" of the smoothing 
kernel should be identical to the resolution of the classical detector. 

There are various possibilities to identify the effective area of a 
two-dimensional centrally peaked function. Probably the most gen- 
eral one is the total dispersion a, defined through 



W{x] 



\x\^&x. 



(19) 



One can readily verify that the classical detector kernel (10) has 
a dispersion a — A/ a/6. Requiring that the resolution of the 
other smoothing kernels are equal to this value we find as reference 
smoothing lengths 



/iref = ^ A 0.408 A 

V6 



/iref ■ 
/iref ■ 



33 
7 



-/\^ 0.598 A 
: A ^ 0.513 A 



(gaussian), (20) 

(Ms spline), (21) 
(M4 spline). (22) 



Figure 1 shows a comparison of the three different smoothing ker- 
nels with their optimal smoothing lengths. 



3 TESTS 

3.1 Comparison of classical and smart detectors 

We tested the accuracy and the performance of our smart detec- 
tors using two toy analytical surface brightness distributions on the 
plane of the sky. The first is a simple Plummer model, characterized 



by the circularly symmetric surface brightness distribution 



(23) 



with L the total luminosity and b sl scale parameter. We took the 
value b = 10 A, implying that the core of the model is well re- 
solved by the detector grid. The second model is an exponential 
disc model, rotated over an angle of 20 degrees. 



7(a 



exp 



\y\ 



(24) 



with hx and hy the scalelength and scaleheight respectively. 
Choosing the disc parameters as hx = 25A and hy = 5A, we 
create a model with a relatively sharp edge and a strong gradient in 
the surface brightness distribution. 

We have used a classical detector and three smart detectors 
based on the smoothing kernels (16), (17) and (18), each of them 
with a total of 101 x 101 grid points. For each model, we ran a set 
of Monte Carlo simulations with the number of photon packages 
varying between = lO'^ and N — 10^, taking into account that 
each detector measures the same set of photon packages (i.e. we 
use the same Monte Carlo realisation for simulations with different 
detectors). 

The different panels of figure 2 show the resulting measured 
surface brightness distribution h{x) for the various detectors for 
the simulations with N — 10^. Looking at this set of images, it 
is immediately clear that the smart detectors manage to qualita- 
tively reproduce the surface brightness distribution accurately. For 
a measure of the accuracy of the different detectors, we consider 
the noise as the difference between the measured surface bright- 
ness and the theoretical surface brightness in each pixel, weighted 
by the expected Poisson noise cr{xij) in each pixel. 



a{xij) 



(25) 



Figure 3 show the noise images corresponding to the Plummer 
model images from Figure 2. It is clear that these images are 
qualitatively very similar; thanks to the appropriate choice of the 
smoothing lengths for the different kernels, there is no correlation 
in the noise on a pixel-by-pixel scale. A quantitative analysis of 
these noise images, however, shows that the level of the noise in 
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Figure 2. The Plummer (top) and exponential disc (bottom) model surface brightness distribution I(xij) and the observed surface brightness distributions 
Is(xij) for the classical and the smart detectors. The images have 101 x 101 pixels each and are based on the simulation with N = 10^ photon packages. 




Figure 3. The noise distributions AI{xij) corresponding to the observed Plummer surface brightness distributions displayed in Figure 2. 



the smart detectors is suppressed. This is most easily demonstrated 
using the total noise level, which we define as 



7^: 



Is^Xij^ I(^Xi 



(26) 



where the sum runs over all pixels. A factor is included in this 
formula to guarantee that the noise is asymptotically independent 
of the total number N of photon packages used in the simulation. 
Figure 4 shows the value of IZ for the different detectors as a func- 
tion of the total number of photon packages in the simulation. This 
figure demonstrates that all smart detectors reduce the noise com- 
pared to the classical detector. The most efficient detector is the 
one based on the M3 spline kernel; for this detector the noise is 
reduced by about 10%. For the M4 spline kernel, the most popu- 
lar finite support kernel in SPH simulations, the noise reduction is 
slightly less efficient (about 8%), whereas for the gaussian kernel 
the noise reduction is about 5%. 



3.2 Origin of the noise reduction 

We have constructed our smart detectors based on two fundamen- 
tal changes applied to the classical detector smoothing kernel (10), 
namely making it circularly symmetric and choosing a smoothly 
decreasing function of radius. We can investigate which of these 



two changes has the most important impact on the noise reduction 
by constructing two new smoothing kernels in which only one of 
these two changes is taken into account. 

On the one hand, the circularly symmetric analogue of the 
classical detector kernel (10) is readily found, 

W{R) = ^[H{R+\h)-H{R-\h)]. (27) 

We find in the usual way /iref = 2 A/ a/3 as reference smoothing 
length. On the other hand, the rectangular version of the M3 spline 
kernel (17) is 



Ms 



with 



M3H 



1 /3 

2 V2 







if < IX < |, 
if| <^< I, 
else. 



(28) 



(29) 



For this kernel we obtain /iref = A/ a/3. 

In Figure 5 we plot the total noise parameter IZ of these two 
new smart detectors in comparison with the classical detector and 
the circularly symmetric M3 spline detector. Not surprisingly, we 
find that also these two new detectors suppress the noise compared 
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Figure 5. The total noise parameter TZ of the Plummer model for the rect- 
angular and circular equivalents of the classical detector and the smart M3 
spline detector. 



to the classical detector, but not as efficiently as the M3 spline de- 
tector. The noise reduction in a detector based on the rectangular 
M3 spline kernel (28) is more effective than the noise reduction in 
the circular analogue (27) of the classical detector. This means that 
applying a proper weight to the photon packages according to their 
distance from the grid points is the most important change if we 
want efficient noise reduction. 



4 DISCUSSION AND CONCLUSIONS 

We have focused on an ill- studied area of Monte Carlo radiative 
transfer simulations where a significant noise reduction can be 
achieved, namely the detection of photon packages and the cor- 
responding construction of the observed surface brightness distri- 
bution. The motivation of this work was the fact that the classical 
detectors used in Monte Carlo simulations, while closely mimick- 
ing a real CCD detector, do not use the full amount of information 
that is available. 

Based on the similarities between the construction of the sur- 
face brightness distribution on a detector in Monte Carlo radiative 
transfer simulations and the calculation of the density in SPH hy- 
drodynamical simulations, we have constructed a set of smart de- 
tectors. These smart detectors improve on two aspects of the clas- 
sical detector: they assign the same weight to all photon packages 
hitting the detector at the same distance from a grid point and they 
give more weight to impacts close to a grid point than to impacts 
at larger distances. We have tested different kinds of smart detec- 
tors, based on different smoothing kernels frequently encountered 
in SPH simulations, namely a gaussian kernel and two spline-based 
kernels with finite support (Hockney & Eastwood 1981; Monaghan 
& Lattanzio 1985). We have shown that these new detectors, while 
preserving the same effective resolution, reduce the noise in the sur- 
face brightness distributions compared to the classical detectors. It 
is demonstrated that the lion's share of this noise reduction is due 
to a proper weighing of the photon packages with the distance be- 
tween the impact location and the grid point. 

The most efficient smart detector is found to be a detec- 
tor based on the M3 spline kernel, for which the noise reduc- 
tion amounts to some 10%. While this might seem a modest im- 
provement, one should take into account that the noise reduction in 
Monte Carlo simulation goes as 1 / A reduction of the noise 



with 10% is hence equivalent to a reduction of the number of pho- 
ton packages with 20%, which is a significant improvement. More- 
over, it should be stressed that this noise reduction basically comes 
for free, since the practical implementation of the smart detectors 
is straightforward and the additional computational cost is com- 
pletely negligible. We hence strongly recommend the use of smart 
detectors in Monte Carlo radiative transfer simulations. 

The link with SPH simulations might stimulate to look for fur- 
ther ways to optimize the estimates of the surface brightness distri- 
bution. One major difference between our current problem and the 
interpolation in SPH simulations is that we use a fixed smoothing 
length, whereas SPH simulations typically use a spatially (and tem- 
porally) varying smoothing length. This way it is possible to take 
full advantage of the particle distribution to resolve local density 
structures. Each particle in an SPH simulation typically has an in- 
dividual smoothing length which is fine-tuned such that each par- 
ticle interacts with a similar number of neighbors. Unfortunately, 
it seems hard to introduce variable smoothing lengths here in our 
current Monte Carlo radiative transfer case. The main difference is 
that the smoothing procedure in SPH is executed when all parti- 
cle positions are known, whereas in Monte Carlo radiative transfer 
the photon packages gradually hit the detector and they must be 
smoothed out onto the detector before it is known where the other 
photon packages will arrive. It is impossible to know a priori how 
many neighbors a photon package will ultimately have and hence 
to adapt the smoothing length on an individual basis. One could po- 
tentially do a test run with a limited number of photon packages to 
obtain a first crude estimate of the expected surface brightness dis- 
tribution and adapt the smoothing lengths accordingly. However, a 
more elegant option seems to apply an adaptive filtering to the sim- 
ulated images after the simulation. Several approaches have been 
developed for this goal such as wavelet-based algorithms (Starck 
& Pierre 1998), adaptive binning techniques (Starck & Pierre 1998; 
Cappellari & Copin 2003) or adaptive kernel smoothing techniques 
(Richter et al. 1991; Huang & Sarazin 1996; Ebeling, White, & 
Rangarajan 2006). 
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